#include "Injection.h"

double GetParticleMass(double rho, double dia)
{
	return (1.0 / 6.0) * M_PI * dia * rho * dia * dia;
}

void GetRRSizeMassRatio(vector<double>& RRsize, vector<double>& RRmassRatio, int RRInterval, double RRMean, double RRSP, double RRmin, double RRmax)
{
	RRsize = vector<double>(RRInterval, 0.0);
	RRmassRatio = vector<double>(RRInterval, 0.0);
	double sizeInterval = (RRmax - RRmin) / (RRInterval - 1.0);

	vector<double> RRsize2(RRInterval, 0.0), Yd(RRInterval, 0.0);

	for (int i = 0; i < RRInterval; i++)
	{
		RRsize[i] = RRmin + sizeInterval * i;

	}
	for (int i = 0; i < RRInterval; i++)
	{
		if (i < RRInterval - 1)
		{
			RRsize2[i] = 0.5 * (RRsize[i] + RRsize[i + 1]);
		}
		else if (i == RRInterval - 1)
		{
			RRsize2[i] = RRsize[i] + 0.5 * (RRsize[i] - RRsize[i - 1]);
		}
		Yd[i] = exp(-1 * pow(RRsize2[i] / RRMean, RRSP));

	}

	vector<double> mass_R2(RRInterval, 0.0), sum_R(RRInterval, 0.0), mass_R(RRInterval, 0.0), mass_R0(RRInterval, 0.0);

	for (int i = 0; i < RRInterval; i++)
	{
		if (i == 0)
		{
			mass_R2[i] = 1.0 - Yd[i];
			sum_R[i] = mass_R2[i];

		}
		else if (i > 0 && i < RRInterval - 1)
		{
			mass_R[i] = Yd[i - 1];
			mass_R0[i] = Yd[i];
			mass_R2[i] = Yd[i - 1] - Yd[i];

		}
		else if (i == RRInterval - 1)
		{
			mass_R2[i] = Yd[i - 1] - Yd[i];
		}
	}
	double totalMass = accumulate(mass_R2.begin(), mass_R2.end(), 0.0);
	for (int i = 0; i < RRInterval; i++)
	{
		RRmassRatio[i] = mass_R2[i] / totalMass;
	}
}

int GetPNumber(shared_ptr<ModelSettings>m1)
{
	/*part 1*/
	auto real_injection_number = (*m1).numberInjectionPosition;
	double injection_interval = ((*m1).endInjectionYPosition - (*m1).startInjectionYPosition) / (real_injection_number + 1);
	vector<double> injection_location(real_injection_number, 0.0);
	for (int i = 0; i < real_injection_number; i++)
	{
		injection_location[i] = (*m1).startInjectionYPosition + injection_interval * (i + 1.0);
	}
	vector<double> injection_area(real_injection_number, injection_interval);

	/*part 2*/
	int size_dist_no;
	if ((*m1).sizeTypes==1)
	{
		size_dist_no = (*m1).PSD_RR_nBins;
	}
	else if ((*m1).sizeTypes == 0)
	{
		size_dist_no = ((*m1).PSD_User_Dia).size();
	}
	int pNumber = real_injection_number * size_dist_no;

	return pNumber;
}

vector<Particle> GetInjectionAll(shared_ptr<CaseSettings> c1, shared_ptr<ModelSettings>m2,
	shared_ptr<Constants>const1)
{
	auto pNumber = GetPNumber(m2);
	vector<Particle> pInjAll(pNumber);

	/*Get injection location and injection area*/
	auto real_injection_number = (*m2).numberInjectionPosition;
	double injection_interval = ((*m2).endInjectionYPosition - (*m2).startInjectionYPosition) / (real_injection_number + 1);
	vector<double> injection_location(real_injection_number, 0.0);
	for (int i = 0; i < real_injection_number; i++)
	{
		injection_location[i] = (*m2).startInjectionYPosition + injection_interval * (i + 1.0);
	}
	vector<double> injection_area(real_injection_number, injection_interval);

	/*Get dicrete particle size and mass percentage ratio based on size*/
	int size_dist_no;
	if ((*m2).sizeTypes == 1)
	{
		size_dist_no = (*m2).PSD_RR_nBins;
	}
	else if ((*m2).sizeTypes == 0)
	{
		size_dist_no = ((*m2).PSD_User_Dia).size();
	}
	vector<double> dia_mass_ratio(pNumber, 0.0), par_inj_area(pNumber, 0.0);
	//vector<int> number_try_index(pNumber, 0);
	vector<double> dia_dist, dia_mass_ratio0;
	if ((*m2).sizeTypes == 1)
	{
		GetRRSizeMassRatio(dia_dist, dia_mass_ratio0, size_dist_no,
			(*m2).PSD_RR_Dmean, (*m2).PSD_RR_SP, (*m2).PSD_RR_Dmin, (*m2).PSD_RR_Dmax);
	}
	else if ((*m2).sizeTypes == 0)
	{
		dia_dist = (*m2).PSD_User_Dia;
		dia_mass_ratio0 = (*m2).PSD_User_DiaRatio;
	}

	/*Loop all particles*/
	for (int i = 0; i < pNumber; i++)
	{
		pInjAll[i].V = vector<double>(2, 0.0);
		pInjAll[i].pos = vector<double>(2, 0.0);

		int index_injection = ceil(double(i + 1.0) / double(size_dist_no));
		int index_in_injection = (i + 1) % size_dist_no;
		if (index_in_injection == 0)
		{
			index_in_injection = size_dist_no;
		}
		index_injection = index_injection - 1;
		index_in_injection = index_in_injection - 1;
		pInjAll[i].diam = dia_dist[index_in_injection];
		dia_mass_ratio[i] = dia_mass_ratio0[index_in_injection];
		pInjAll[i].pos[1] = injection_location[index_injection];
		par_inj_area[i] = injection_area[index_injection];
		pInjAll[i].ID = i;

		pInjAll[i].time = 0.0;
		pInjAll[i].V[0] = (*m2).pVx0;
		pInjAll[i].V[1] = (*m2).pVy0;
		pInjAll[i].pos[0] = (*c1).coordinateVertex0[0];
		pInjAll[i].temp = (*m2).pT0;
		pInjAll[i].rho = (*m2).pRho;
		pInjAll[i].mass = GetParticleMass(pInjAll[i].rho, pInjAll[i].diam);
		/*Get SF*/
		double total_par_inj = (*m2).totalMassFlowRate; // unit, kg/s
		double mass_par = total_par_inj * dia_mass_ratio[i];
		pInjAll[i].SF = mass_par / pInjAll[i].mass;
	}

	return pInjAll;
}

void debugInjection(vector<Particle> pAll)
{
	string injectionFileCheck = "debugInjection0.txt";
	ofstream write_injection(injectionFileCheck);
	assert(write_injection.is_open());

	int length = pAll.size();

	for (int i = 0; i < length; i++)
	{
		write_injection.precision(10);
		if (i == 0)
		{
			write_injection << "ID" << "," << "V[0]" << "," << "V[1]" << "," << "pos[0]" << "," << "pos[1]"
				<< "," << "diam" << "," << "rho" << "," << "mass" << "," << "temp" << "," << "time"
				<< "," << "SF" << "\n";
		}
		write_injection << pAll[i].ID<<","<< pAll[i].V[0] << "," << pAll[i].V[1] << "," << pAll[i].pos[0] << "," << pAll[i].pos[1]
			<< "," << pAll[i].diam << "," << pAll[i].rho << "," << pAll[i].mass << "," << pAll[i].temp << "," << pAll[i].time 
			<< "," << pAll[i].SF << "\n";
	}
	write_injection.close();
}